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Summary 

Mathematical  formulae  and  a  Fortran  program  listing  are 
given  for  calculating  the  response  and  sound  radiation  due  to 
prescribed  time-harmonic  excitation.  The  program  is  useful 
for  bench-mark  studies.  Numerical  results  illustrate  the 
applicability  of  the  various  ■''DAA'’'”  forms  of  the  exterior 
fluid  loadincr  term.  j  „ 
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1.  INTRODUCTION 


The  submerged  spherical  shell  is  the  only  finite  elastic  system 
for  which  closed-form  expressions  are  available  for  its  vibration 
and  sound  radiation.  Its  simplicity  makes  it  a  useful  bench-mark 
problem  against  which  purely  numerical  formulations  may  be 
tested.  Much  of  the  theoretical  work  on  the  dynamics  and 
acoustics  of  the  spherical  shell  is  contained  in  the  standard 
text  C13  on  fluid-structure  interaction,  authored  by  Junger  and 
Feit.  The  vibration  of  a  point  force  excited  spherical  shell  in 
air  and  in  water  has  been  studied  by  Mong  &  Hayek  C23  who  found 
good  agreement  between  the  resonant  frequencies  determined  by 
theory  and  experiment.  James  C33  has  presented  target  strength 
spectra  and  intensity  vector  plots  which  show  the  considerable 
distortion  of  the  sound  field  in  the  vicinity  of  a  spherical 
shell,  when  it  is  excited  at  a  resonant  frequency.  Ansley  & 
Skelton  C43  have  studied  the  scattering  of  sound  by  multilayered 
spherical  shells  using  the  exact  theory  of  linear  elasticity: 
some  of  their  numerical  work  (not  published)  demonstrates  that 
the  shell  theory  approximation  is  valid  for  radiated  sound 
calculations  over  a  wide  range  of  frequency  and  thickness  to 
radius  ratios. 

Recent  theoretical  research  in  f luid-structure  interaction 
problems  has  developed  approximations  which  avoid  the  lengthy 
computations  required  by  the  Helmholtz  integral  equation.  These 
theories,  which  are  collectively  known  by  the  name  'doubly 
asymptotic  approximations  (DAA's)',  converge  correctly  at  low  and 
high  frequencies  to  the  'added  mass'  and  'plane  wave'  approx¬ 
imations,  respectively.  They  have  been  obtained  by  Geers  C53 
and  summarized  by  Geers  &  Felippa  £63. 

The  vibration  and  sound  radiation  of  a  spherical  shell  is 
the  standard  bench-mark  problem  for  testing  the  performance  of 
the  various  DAA's,  both  analytically  and  numerically.  Geers  & 
Felippa  £63  have  given  closed-form  expressions  for  the  DAA's  of 
the  submerged  spherical  shell,  and  have  presented  numerical 
results  which  show  that  the  second  order  approximations  are 
adequate,  except  at  frequencies  close  to  resonances.  Their 
findings  have  been  reinforced  by  the  numerical  calculations  of 
Huang  &  Wang  £73. 

The  numerical  studies  of  Geers  &  Felippa  are  limited  to 
frequency  response  plots  of  the  shell's  velocity  and  surface 
pressure  in  the  individual  terms  (n=0,3)  of  their  Legendre  series 
expansions.  The  numerical  comparisons  of  Huang  &  Wang  are  based 
on  polar  plots  of  the  far-field  sound  radiation,  at  selected 
frequencies:  their  comparisons  could  be  misleading  close  to  a 
resonant  frequency  because  small  differences  in  the  resonant 
frequencies,  which  in  themselves  are  of  no  practical  importance, 
could  cause  unduly  large  differences  between  the  sound  pressure 
levels  at  a  selected  frequency.  Thus,  in  addition  to  providing 
a  computer  program  which  can  be  used  for  bench-mark  studies,  it 
is  the  purpose  of  this  Memorandum  to  extend  the  range  of  numer¬ 
ical  tests  of  the  applicability  of  the  DAA's  to  the  spherical 
shell  problem. 


!.  RESUME  OF  MATHEMATICS 

The  radial  and  tangential  displacements  of  the  shell's 
surface  are  expressed  as  a  series  of  Legendre  polynomials  of 
argument  vt*cos6,  viz. 


U(0)  =  Z  (l-p2)1/2UnPn(u) 
n=0 


(1) 


W(0)  *  Z  W  P  (y) 
n  n 

n=0 
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whose  amplitudes  are  obtained  from  the  matrix  relation 
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where  for  coupled  membrane  and  bending,  theory  in  which  effects 
due  to  rotary  inertia  and  transverse  shear  are  neglected 
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a22 

=  E.P2x  (cr+x  -1)+2E.  (l+a)-u)2p  h 
inn  1  s 

+p2C2u,hn(k2a)/hn(k2a)  "  plclwVkla  >  'V  kla } 


In  the  above,  h  is  the  shell's  thickness  and  a  its  mean  radius; 

2  2 

E^=Eh/(l-cr  )a  where  E  is  Young's  modulus  and  a  is  Poisson's 
2  2  2 

ratio;  0  =h  /12a  and  x  =n(n+l);  p  is  the  density  of  the 
shell’s  material,  and  p1  and  p2  are  the  densities  of  the 
interior  and  exterior  fluids  whose  sound  velocities  are  c^  and 
c2;  k^  and  k2  are  the  acoustic  wavenumbers,  w/c^  and  w/c2; 
and  hn(=jn+iyn)  are  spherical  Bessel  functions;  Pn  are  Legendre 

polynomials  of  degree  n;  the  prime  on  the  various  quantities 
denotes  differentiation  with  respect  to  their  arguments.  The 
time-harmonic  factor  exp(-iwt)  is  omitted  throughout.  The 
geometry  of  this  problem  is  shown  in  Figure  1,  where  in  a  system 
of  spherical  coordinates  <R,0,$)  the  field  quantities  of  interest 
are  independent  of  the  circumferential  angle,  0. 

The  acoustic  pressure  in  the  exterior  fluid  is  given  by  the 
expression 
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p(r,0)  =  p_c_w  Z  W  P  (cos0)h  (k~R) /h  (k_a) 
2  2  nn  n  2  n  2 

n  =  0 


(5) 


and  the  far-field  pressure  is  obtained  from  this  equation  by 
replacing  the  spherical  Hankel  function  by  its  value  for  a  large 
argument,  giving 


pf(R,0)  =  p2c2wCexp( ik2R) /k2R3  x 


CO 

I  W  P  ( cos0)  expC-i  ( n+1 )  ir/23/h'  (  k_a ) 
n  n  n  2 

n=0 
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3.  EXCITATION 

In  equation  (3)  the  coefficients  Fn  of  the  Legendre 

series  expansion  of  the  prescribed  radial  stress  distribution, 
F(0),  are  defined  by  the  transform  pair 


F( 6 )  =  Z  F  P  <u>  (7) 

n  n 

n  =  0 


F  =  ( n+i / 2 ) ft  F(p)P  (u)du  (8) 

n  - 1  n 

For  the  particular  case  of  a  radial  point  force  FQ  excitation 
at  0=0, 


F(  6 )  =  FQS(0) /2ira2sin0 

F  =  (Frt/2ira2)  (n+1/2 ) 
n  0 


(9) 


Additionally,  if  the  external  stress  is  a  constant  FQ  over  the 
axisymmetric  region  defined  by  0=0  to  a  and  $=0  to  2ir,  then 

Fn  =  (F0/2)CPn_1(cosoi)-Pn+1<cos(x)3  (10) 

in  which  This  is  the  excitation  used  in  reference  C7D. 


4.  EXTERIOR  FLUID  LOADING 

The  exact  relation  between  the  modal  acoustic  pressure  and 
the  radial  displacement  at  the  shell's  surface  is 

p„  =  pcwth  (ka)/h  <ka)DW„  (11) 

n  n  n  n 

where  p  and  c  are  the  density  and  sound  velocity  of  the  exterior 
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fluid,  and  k=w/c.  In  the  shell's  equations  (3)  and  (4)  the 
effect  of  the  surrounding  fluid  is  contained  in  the  coefficient 
a22  of  the  matrix  relation,  equation  (3), 


If  S(n,w)  is  defined  as  the  modal  receptance  of  the 
displacement  and  pressure  at  the  shell's  surface,  then  its 
various  analytical  approximations  can  be  found  in  Geers  &  Felippa 
C63  and  Huang  &  Hang  C73  as: 


(0)  exact 

# 

S(n,w)  =  pcwh  (ka)/h  (ka) 
n  n 


(1)  plane  wave  approximation,  PWA 
S(n,w)  =  -iwpc 


(2)  added  mass  approximation,  AMA 

2 

S(n,w)  =  -paw  /(n+1) 


(3) 


doubly  asymptotic  approximation,  OAA1 
2 

S(n,w)  =  paw  /Cika-(n+l)3 


( 12) 


(4) 


doubly  asymptotic  approximation,  DAA2(M) 

S(n,w)  =  paw2Cika- ( n+1 )  3/C  ( n+1 ) 2-ika( n+1 )  -k2a2D 


(5) 


doubly  asymptotic  approximation,  DAA2(C) 
S(n,w)  =  paw2Cika-n3/Cn(n+l)-ika(n+l)-k2a23 


In  Geers  &  Felippa  C63,  DAA2(M)  is  described  as  a  mode  derived 
DAA2,  and  DAA2(C)  is  described  as  a  curvature  augmented  DAA2. 

Thus,  if  the  terms  of  equations  (12)  are  used  in  turn  in  equations 
(4),  numerical  plots  of  vibration  and  sound  radiation  can  be  used 
as  a  test  of  the  range  of  validity  of  the  various  fluid  loading 
approximations  of  this  problem. 


5.  NUMERICAL  RESULTS 
( a )  General 


A  Fortran  program  has  been  written  to  calculate  radial 
response  and  the  far-field  sound  radiation  of  the  point  force 
excited  shell.  A  listing,  together  with  the  computer  program 
specifications  and  test  examples,  is  given  in  the  Appendices. 

The  computational  procedure  for  the  far-field  sound  radiation  is 
straightforward,  because  the  Hankel  function  in  the  denominator 
ensures  rapid  convergence  of  the  series  in  equation  (6).  The 
series  for  the  displacement  in  equation  (2)  converges  slowly,  so 
in  this  case  the  fluid  loading  terms  for  n>>ka  are  best  computed 
from  their  asymptotic  form 

-P2c2wk2a/ (n+1 )  -  p^c^uk^a/n  (13) 
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with  negligible  loss  of  accuracy. 


In  Figures  2  and  3  which  show  spectra  of  the  shell's  dis¬ 
placement  and  far-field  pressure  at  0=0  for  the  exact  and  the 
various  exterior  fluid  approximations,  the  following  constants  in 
SI  units  were  used: 

Steel  shell:  E=19.5xl010  o= 0.29  p  =7700  a=1.0  h=0.01 

s 

Exterior  water:  p=1000  c=1500  Interior  vacuum:  p=0  c=0 

Shell  damping  was  included  in  the  computations  by  setting  Young's 
modulus  to  the  complex  value  E=E(i-in>,  where  the  hysteretic 
loss-factor  p  was  chosen  as  0.01.  In  the  plots,  the  shell's 
displacement  is  expressed  in  dB  ref.  1  micrometre,  and  the 
pressure  is  in  dB  ref.  1  micropascal  at  lm.  The  excitation  is  a 
unit  force  applied  at  0=0. 

(b)  Shell  Displacement 

The  displacement  for  the  EXACT  exterior  fluid  loading  is 
shown  in  Figure  2A.  At  20Hz  the  anti -resonance  is  caused  by 
interaction  between  the  n=l  rigid  body  mode  and  the  forced 
response  of  the  other  modes  excited  well  below  their  resonant 
frequencies,  which  are  defined  herein  as  the  frequencies  at  which 
the  modal  amplitudes  are  a  maximum.  Above  250Hz  the  shell's 
resonances  dominate  in  the  response,  the  resonance  spacing  being 
a  minimum  near  to  500Hz.  Above  this  frequency  the  shell  becomes 
increasingly  stiff  as  its  mean  displacement  limits  to  that  of  a 
water  loaded  infinite  plate. 

Figure  2B  shows  the  displacement  for  the  PWA  fluid  loading: 
this  approximation  is  clearly  inadequate  at  all  frequencies. 

The  AMA  plot  of  Figure  2C  and  the  EXACT  plot  are  in  good 
agreement,  except  for  the  n=2-4  modes  whose  AMA  amplitudes  are 
too  high  and  too  far  shifted  to  higher  frequencies.  The  DAA1 
plot  of  Figure  2D  agrees  well  with  the  EXACT  plot  up  to  250Hz, 
when  the  excessive  radiation  damping  inhibits  the  formation  of 
resonances.  Figures  2E  &  2F  show  that  the  DAA2(M)  and  DAA2(C> 
plots  agree  with  the  EXACT  plot  in  all  aspects,  except  that  the 
resonant  peaks  are  underestimated  by  a  maximum  of  some  4  dB. 

( c )  Far-field  Pressure 


The  far-field  pressure  for  EXACT  exterior  fluid  loading  is 
shown  in  Figure  3A.  At  frequencies  up  to  200Hz  the  pressure  is 
close  to  that  obtained  from  a  point  force  on  an  infinite  plate. 
Between  200  and  500Hz  the  spectrum  is  dominated  by  a  sharp 
resonant  structure.  Above  500Hz,  the  pressure  is  again  close  to 
the  pressure  of  an  infinite  plate.  In  this  frequency  regime  the 
radiation  resistance  of  the  modes  at  resonance  is  too  small  for 
them  to  contribute  to  the  far-field  pressure,  which  is  due  to  the 
radiation  from  modes  excited  above  their  resonant  frequencies. 

Figure  3B  shows  the  sound  radiation  for  the  PWA  fluid 
loading:  this  approximation  is  poor  except  at  frequencies  above 
900Hz ,  where  the  radiation  of  the  EXACT  plot  is  due  to  modes 
radiating  with  pc  efficiency.  The  AMA  plot  of  Figure  3C  has 


much  the  same  shape  as  the  EXACT  plot,  but  the  resonant  peaks  are 
too  high  and  too  far  shifted  to  higher  frequencies:  also,  the 
spectrum  is  3dB  too  low  above  500Hz.  Except  at  very  low 
frequencies  the  DAA1  plot  of  Figure  3D  is  a  poor  approximation. 
The  levels  of  the  n=3-5  modes  of  the  DAA2(M)  plot  of  Figure  3E 
are  some  lOdB  too  low,  but  their  frequency  location  is  much  the 
same  as  in  the  EXACT  plot:  both  low  and  high  frequency  regimes 
provide  good  agreement  with  those  of  the  EXACT  plot.  The  DAA2(C) 
plot  of  Figure  3F  agrees  less  well  with  the  EXACT  plot  than  does 
the  DAA2(M)  plot. 


6 .y CONCLUDING  REMARKS 

The  principal  features  of  the  work  contained  herein  are 

these: 

(1)  A  Fortran  computer  program  which  can  be  used  for 
bench-mark  studies  of  purely  numerical  formulations 
of  the  vibration  and  sound  radiation  of  submerged 
elastic  structures. 

(2)  Numerical  results  of  vibration  and  sound  radiation 
which  augment  those  obtained  from  other  published 
studies,  and  reinforce  the  view  that  better  fluid 
loading  approximations  are  required  in  the  inter¬ 
mediate  frequency  regime  where  resonances  occur,  ■..e*" 

From  a  practical  point  of  view  it  may  be  thought  that  AMA  and 
DAA2(M)  give  acceptable  results  for  both  vibration  and  sound 
radiation,  especially  if  the  application  is  more  suited  to 
phenomenological  rather  than  one-to-one  modelling,  which  is 
usually  the  case  in  noise  control  engineering. 


I 


J.H.  James  (PSO) 
November  1985 
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PROGRAM  SPECIFICATIONS 


(a)  Data  Preparation 


Item  01;  6  real  numbers,  describing  the  shell. 


ES 

VS 

RHOS 

AS 

HS 

ETAS 


Young's  modulus 
Poisson's  ratio 
density 
radius 
thickness 

hysteretic  loss  factor 


Item  02;  4  real  numbers,  describing  the  fluids. 


RHOl 

Cl 

RH02 

C2 


density  interior  fluid 
sound  velocity  interior  fluid 
density  exterior  fluid 
sound  velocity  exterior  fluid 


Item  03:  1  integer,  2  reals  and  3  integers. 


NF 

FS 

FM 

NHS 

NHM 

NEX 


number  of  frequencies  in  linear  frequency  sweep 
start  frequency 
final  frequency 

number  of  harmonics  summed  at  first  frequency 
number  of  harmonics  summed  at  final  frequency 
number  of  extra  harmonics  summed  for  shell  response 


Item  04:  2  integers,  1  real  and  1  integer. 


INTI 

IND2 

THDEG 

IPLOT 


=0,1, 2, 3  for  computation  of  far-field  pressure, 

shell  displacement,  velocity,  acceleration 

=0,1, 2, 3, 4, 5  for  EXACT, PWA, AMA,DAA1 ,DAA2(M) ,DAA2(C) 

computation  of  exterior  fluid  loading 

angle  in  degrees  where  response  is  required,  the 

drive  point  being  at  THDEG=0.0 

store  on  disk  in  block  IPLOT.  Machine  dependent. 


(b)  Notes  on  Program 

(1)  In  Item  02,  the  value  of  the  sound  velocity  can  be  set  to 
zero  if  the  density  is  given  as  zero. 

(2)  In  Item  03,  the  start  frequency  must  be  greater  than  zero,  a 
value  of  about  0.01c /a  being  suitable.  The  values  of  NHS  and 
NHM  should  be  selected  by  trial  and  error:  values  too  large 
will  cause  overflow.  NHM  should  exceed,  by  at  least  6,  say, 
the  number  of  resonances  in  the  displacement  response,  and 

it  should  also  be  greater  than  5ka.  The  number  of  harmonics 
summed  at  intermediate  frequencies  is  obtained  by  linear 
interpolation. 


(3)  The  integer  NEX  in  Item  03  has  no  effect  if  the  far-field 


pressure  is  required.  It  is  needed  for  the  shell's  response 
because  of  the  slow  convergence  of  the  series.  A  value  such 
that  (NEX+NHM)=90  appears  to  be  suitable.  If  (NEX+NHM)>98 
the  array  dimensions  in  SUBROUTINE  VIBRAD  must  be  increased. 
The  number  of  frequencies  NF  must  not  exeed  501  unless  the 
array  dimensions  in  the  MAIN  program  are  increased. 

(4)  If  IND1=0  in  Item  04  the  far-field  pressure  is  given  in  dB  ref 
1  micropascal  at  lm  for  unit  excitation  force,  SI  units  being 
assumed.  If  IND1=1,2,3  the  shell  displacement,  velocity  or 
acceleration  is  given  in  dB  ref.  1  micro-unit  for  unit  excit¬ 
ation  force. 

(5)  The  integer  IND2  in  Item  04  is  for  the  approximation  (if  any) 
required  for  the  exterior  fluid  loading  term.  IND2=0  gives 
the  exact  solution,  and  IND2=1  to  5  gives  the  approximations 
of  Section  4.  The  interior  fluid  loading  is  always  exact. 


c)  Test 

Example 

for  Sound  Radiation 

9.5E10 

0.29  7700.0  1.0 

0.01 

0.01 

i.O  0.0 
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c  APPENDIX  B 

C . 

C  COPYRIGHT  CONTROLLER  HER  MAJESTY'S  STATIONARY  OFFICE 
C  LONDON,  1986. 

C . 

c 
c 

PROGRAM  S HR AD 2 
C 

C  JH  JAMES.  ARE ( TEDDINGTON ) . 

C  VIBRATION  AND  RADIATION  OF  FLUID  FILLED  SPHERICAL  SHELL 
C 

COMPLEX  VALUE 

DIMENSION  DB( 501 ) ,PHASE( 501 ) ,HZ( 501 ) 

COMMON / SHI / ES , VS , RHOS , AS , HS , ETAS 
COMMON /SH2/RH01 ,C1 ,RH02 ,C2 ,FREQ 
COMMON / SH3 / IDV1 , IDV2 , NHARM , NEX 
C 

600  FORMAT( II, IX, '  FREQ','  DB',2X,'  PHASE') 

602  FORMAT ( IX) 

604  FORMAT ( IX ,F9 . 1 ,F7 . 1 , 2X,F7 . 1 ) 

C 

IDV1=5 

IDV2=7 

READ( IDV1 , * )ES , VS , RHOS , AS ,HS ,ETAS 
READ( IDV1 , * )RH01 ,C1 ,RH02 ,C2 
READ ( IDV1 ,*) NF , FS , FM , NHS , NHM , NEX 
READ( IDV1 , * ) IND1 , IND2 .THDEG , I PLOT 
WRITE( IDV2 ,600 ) 

C 

C  LOOP  FOR  EACH  FREQUENCY. 

DO  30  1=1, NF 
NFREQ=I 

CON=FLOAT( NFREQ-1 ) /FLOAT( NF-1 ) 

FREQ=FS+ ( FM-FS ) *CON 
NHARM=NHS+ ( NHM-NHS ) *CON 
R=1 . 0 

CALL  VIBRAD(R, THDEG, IND1,IND2, VALUE) 

IF(MOD(NFREQ,10) .EQ.1)WRITE( IDV2,602) 

XR=REAL( VALUE) 

XI =AIMAG( VALUE) 

DBVAL=10.0*ALOG10(XR**2+XI**2+1.0E-30) +120.0 
PHVAL=ATAN2 ( XI ,XR ) *57 . 296 
WRITE( IDV2 ,604 ) FREQ , DBVAL , PHVAL 
HZ  <  NFREQ ) =FREQ 
DB ( NFREQ ) = DBVAL 
PHASE ( NFREQ ) = PHVAL 
30  CONTINUE 

IF ( IPLOT . LE . 0 ) STOP 
C 

C  STORE  ON  DISK  FILE  FOR  PLOTTING.  SYSTEM  DEPENDENT. 

CALL  ASSIGN ( 1, 'DK:FTN1 .DAT' ,0, 'OLD' ) 

DEFINE  FI LEI ( 8, 4 096, U, I VAR) 

WRITE ( 1 ' IPLOT )NF, (HZ( I ) ,DB( I ) , PHASE ( I ) , 1*1 ,NF) 

STOP 

END 

C 

C 
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SUBROUTINE  VIBRAD(R,THDEG,IND1 , IND2, VALUE) 

C 

C  FAR -FIELD  PRESSURE,  AT  DISTANCE  R,  OR  SHELL  RESPONSE. 

C 

REAL  K1,K2 

COMPLEX  VALUE , FLUID , WN , CXI , CX2 , CX3 , CXZERO , CXI 
DIMENSION  AJ1 ( 101 ) ,AY1 (101) ,AJ1D( 101 ) , AY1D( 101 ) 

DIMENSION  AJ2 ( 101 ) ,AY2 ( 101 ) ,AJ2D( 101 ) ,AY2D( 101 ) 

DIMENSION  PN( 101 ) ,PND( 101 ) 

COMMON/ SHI / ES , VS , RHOS , AS , HS , ETAS 
COMMON / SH2 / RHO 1 , C 1 , RHO  2 , C  2 , FREQ 
COMMON / SH3 / IDV1 , IDV2 , NHARM , NEX 
C 

NDIM=101 
PI=3. 141592654 
CXZERO* ( 0 . 0 , 0 . 0 ) 

CXI=(0.0,1.0) 

W=2.0*PI*FREQ 

IF ( RHOl . EQ . 0 . 0 ) GO  TO  10 

K1=W/C1 

ARG1=K1*AS 

CALL  SPJY ( AJ1 ,AY1 ,AJ1D,AY1D,NDIM,ARG1 ) 

10  IF ( RH02 . EQ .0.0) GO  TO  12 
K2=W/C2 
ARG2=K2*AS 

CALL  SPJY ( AJ2 ,AY2 , AJ2D,AY2D,NDIM,ARG2 ) 

ARG3=K2*R 

12  ARG4=C0S( THDEG*PI / 180 . 0  > 

CALL  LEGEND (PN, PND, NDIM,ARG4, NHARM) 

C 

C  SUM  RADIAL  SHELL  RESPONSE  OR  FAR  FIELD  PRESSURE. 

VALUE=CXZERO 

DO  20  NSUM=0, NHARM 

N=NSUM 

FLUID=CXZERO 

IF( RH02 . EQ .0.0 )GO  TO  22 

CX1=CEXP( CXI* (ARG3 -FLOAT ( N+l ) *PI / 2 . 0 ) ) /ARG3 
CX2=CMPLX( AJ2D( N+l ) ,AY2D(N+1) ) 

CX3=CMPLX( AJ2 ( N+l ) ,AY2(N+1) ) 

FLUID*RH02*C2*W* ( CX3/CX2 ) 

IF ( IND2 . NE . 0 ) CALL  APPR OX ( FLU ID,N,RH02,C2,W,AS,IND2) 

22  IF(RH01 .NE. 0 . 0 )FLUID*FLUID-RH01*C1*W* ( AJ1 (N+l ) /AJ1D<  N+l ) ) 
CALL  SHELL (N,W, FLUID, WN) 

WN=WN*FLOAT<  2*N+1 > / ( 4 . 0*PI*AS**2 ) 

IF( IND1 .NE. 0 ) VALUE* VALUE+WN*PN( N+l ) 

IF( IND1 . EQ . 0 ) VALUE=VALUE+RH02*C2*W*WN*PN ( N+l ) *CX1 /CX2 
20  CONTINUE 

IF( IND1 . EQ . 0 ) RETURN 
IF ( NEX .  LE .  0 )  GO  TO  50 
C 

C  SUM  ADDITIONAL  HARMONICS  FOR  SHELL  RESPONSE. 

N1=NHARM+1 

N2=N1+NEX-1 

CALL  LEGEND (PN, PND, NDIM,ARG4,N2) 

DO  30  NSUM*N1,N2 

N=NSUM 

FLUID=CXZERO 

IF ( RH02 . EQ . 0 . 0 ) GO  TO  32 


non  no  o  oooo  no 


FLUID=RH02*C2*W*< -ARG2 /FLOAT (N+l ) ) 

IF( IND2.NE. 0)CALL  APPROX <  FLUID, N ,RH02 ,C2 ,W, AS , IND2 ) 
32  IF ( RHO 1 . NE . 0 . 0 ) FLUID=FLUID-RH01*C1*W* ( ARG1 / FLOAT <  N  >  ) 
CALL  SHELL (N,W, FLUID, WN) 

WN=WN*FLOAT( 2*N+1 ) / ( 4 . 0*PI*AS**2 ) 

VALUE=VALUE+WN*PN ( N+l ) 

30  CONTINUE 

50  IF( IND1 .EQ . 2 )VALUE=-CXI*W* VALUE 
IF( IND1.EQ. 3 ) VALUE* -W*  *  2  * VALUE 
RETURN 
END 


SUBROUTINE  APPROX < FLUID ,N,RHO,C,W,A,IND) 

APPROX.  FLUID  LOADING  FOR  SPHERICAL  SHELL. 
IND=l/2/3/4/5  FOR  PWA/AMA/DAA1 /DAA2 < M) /DAA2 ( C ) 

COMPLEX  FLUID, CXI 


CXI=(0. 0,1.0) 

FN=FLOAT(N) 

FNI =FLOAT ( N+l ) 

C0N=RH0*A*W**2 

AK=W*A/C 

IF( IND.EQ. 1 )FLUID=-CXI*W*RHO*C 
IF( IND.EQ. 2 ) FLUID* -CON/ FN1 
IF ( IND.EQ . 3 ) FLUID=CON/ ( CXI*AK~FN1 ) 
IF ( IND . EQ . 4 ) FLUID=C0N* ( CXI+AK-FN1 ) / 
1  (FN1**2-CXIAAK*FN1-AK**2) 

IF  < IND . EQ . 5 ) FLUID=C0N* ( CXI AAK-FN ) / 

1  (FN*FN1-CXI*AK*FN1-AK**2) 

RETURN 

END 


SUBROUTINE  SHELL ( N,W, FLUID, WN) 

SHELL  DISPLACEMENT  FOR  UNIT  RHS . 

COMPLEX  FLUID ,WN,E1,A11,A12,A21,A22, DET 
COMMON / SHI / ES , VS , RHOS , AS , HS , ETAS 

W2=W**2 

El =ES*HS*CMPLX (1.0, -ETAS ) / ( 1 . 0 -VS*VS  > 
E1=E1/AS**2 

BErA2=HS**2/(12.0AAS*A2) 

DN=FLOAT(NA(N+l) ) 

X1=VS+DN-1 . 0 
X2=l . 0+VS 

A1 1 = ( 1 . 0+BETA2 ) *X1 -W2ARH0S AHS /El 

A12=BETA2*X1+X2 

A21*DNAA12 

A22*BETA2ADNAXl+2 . 0*X2-W2*RHOS*HS/E1+FLUID/E1 

DET=A11AA22-A12*A21 

WN=(A11/DET)/E1 

RETURN 
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non  on  o  ono  no 


UNLIMITED 


END 


SUBROUTINE  SPJY(AJ,AY,AJD,AYD,NDIM,X> 

EVALUATES  THE  SPHERICAL  BESSEL  FUNCTIONS  AND  DERIVATIVES. 

DIMENSION  AJ(NDIM) ,AY(NDIM) ,AJD(NDIM) ,AYD(NDIM) 

COMMON / SH3 / IDV1 , IDV2 , NHARM , NEX 


N=2*NHARM+1 
IF(N.GE.NDIM)GO  TO  70 
AJ  (N) =0. 0 
AJ ( N-l ) =1 . 0E-20 
IL=N-2 

DO  40  1=1, IL 
IP=N-1-I 

40  AJ ( IP) = FLOAT ( 2*IP+1 ) *AJ ( IP+1 ) /X-AJ ( IP+2 ) 
T=(SIN(X)/X)/AJ(1) 

DO  50  1=1, N 
50  AJ(I)=AJ(I)*T 
AY ( 1 ) =-C0S (X) /X 
AY ( 2 ) =-COS (X) / (X**2 ) -SIN(X) /X 
N=NHARM 
DO  60  1=1, N 

60  AY( 1+2 ) =FLOAT ( 2*1+1 )*AY< 1+1 ) /X-AY( I ) 

A  JD ( 1 ) = - A J ( 2 ) 

AYD( 1 ) =-AY ( 2 ) 

DO  80  1=1, N 

AJD( 1+1 ) = ( FLOAT ( I ) *AJ ( I ) -FLOAT ( 1+1 ) *AJ( 1+2 ) ) / FLOAT ( 2*1+1 ) 
80  AYD( 1+1 ) =( FLOAT! I ) *AY< I ) -FLOAT ( 1+1 ) *AY( 1+2 ) ) / FLOAT ( 2*1+1 ) 
RETURN 

70  WRITE( IDV2 ,600 ) 

600  FORMAT!//, IX, 'REDUCE  FREQUENCY  OR  INCREASE  ARRAY  SIZE') 
STOP 
END 


SUBROUTINE  LEGEND (PN, PND, NDIM,X,N> 

LEGENDRE  POLYNOMIAL  AND  DERIVATIVE.  MAXIMUM  ORDER  N. 

DIMENSION  PN(NDIM) ,PND(NDIM) 

PN ( 1 ) =1 . 0 
PN ( 2 ) =X 
NM1=N-1 
DO  10  1=1, NM1 

PN( 1+2 ) = (X*FLOAT( 2*1+1 ) *PN( 1+1 ) -FLOAT ( I ) *PN< I ) ) /FLOAT! 1+1 ) 
10  CONTINUE 
PND ( 1 ) =  0 . 0 
DO  20  1=1, N 

20  PND! 1+1 ) =X*PND( I ) +FL0AT! I ) *PN ( I ) 

RETURN 


